# 布鲁氏菌病发病率地图绘制 ----------------------------------------------------------
## 代码结构：
### 1. 环境准备 - 加载包和参数
### 2. 数据准备 - 读取地理数据和发病率数据
### 3. 数据处理 - 计算指标并匹配地理数据
### 4. 可视化设计 - 构建地图可视化框架
### 5. 成果输出 - 保存高质量图片
### 6. 质量验证 - 数据完整性检查

# 1. 环境准备 ------------------------------------------------------------------
## 加载必要的包
library(sf)       # 地理数据处理
library(dplyr)    # 数据清洗转换
library(tmap)     # 专题地图绘制

## 设置地理数据处理选项
tmap_options(check.and.fix = TRUE)  # 自动修复无效地理多边形


# 2. 数据准备 ------------------------------------------------------------------
## 读取新疆地州级行政区划地理数据（GeoJSON格式）
map_path <- "./data/origin/新疆维吾尔自治区_地州.json"
geo_data <- st_read(map_path) %>% 
  st_make_valid() %>%   # 修复可能存在的无效多边形
  filter(!st_is_empty(geometry))  # 移除非空地理单元

## 加载发病率计算数据（来自外部脚本）
source("./code/2.3 人口数估发病率.R")  # 加载 rate_year_city 数据集


# 3. 数据处理 ------------------------------------------------------------------
## 计算2005-2023年累积发病率指标
cumulative_incidence <- rate_year_city %>%
  filter(between(year, 2005, 2023)) %>%  # 筛选目标年份
  group_by(地区, 地区编码) %>%          # 按地州分组
  summarise(
    total_count = sum(count, na.rm = TRUE),   # 总病例数（处理缺失值）
    total_population = mean(value, na.rm = TRUE),  # 总暴露人口（人年数）
    .groups = "drop"
  ) %>% 
  mutate(
    incidence_rate = (total_count / total_population) * 10^5  # 计算每10万发病率
  )
cumulative_incidence$incidence_rate %>% summary()

cumulative_incidence$incidence_rate %>% unique()

## 匹配地理数据与发病率数据
geo_incidence <- geo_data %>% 
  left_join(cumulative_incidence, by = c("name" = "地区")) %>%  # 按名称匹配
  mutate(
    incidence_rate = replace_na(incidence_rate, 0)  # 缺失值填充为0发病率
  ) %>% left_join(nameCE_city, by = c("地区编码" = "city_code"))

# 4.1 可视化设计（发病率） -----------------------------------------------------------------
## 构建基础地图框架
epidemic_map <- tm_shape(geo_incidence, bbox = st_bbox(geo_incidence)) +  # 设置绘图范围
  
  ## 核心制图元素：发病率分级填色
  tm_polygons(
    col = "incidence_rate",        # 映射到发病率字段
    title = "Cumulative Incidence\n(per 100,000, 2005-2023)",
    palette = "YlOrRd",  # 使用WHO推荐的色阶
    style = "fixed",     # 使用固定间隔分级
    breaks = c(0, 100, 200, 300, 400, 500, 600, 700, 800,900), # 根据数据范围设置分组
    border.col = "gray40", # 边界颜色
    lwd = 0.6,           # 边界线宽
    alpha = 0.8          # 透明度
  ) +
 
  ## 地图布局设计
  tm_layout(
    # main.title = "新疆地州级布鲁氏菌病累积发病率时空分布",  # 主标题
    # main.title.size = 1.1,        # 标题字号（相对比例）
    # legend.outside = TRUE,         # 图例外置
    # legend.outside.position = "right",  # 图例右侧外置
    legend.position =  c("left", "top"),
    legend.format = list(
      fun = function(x) format(x, digits = 2)  # 图例数值保留2位小数
    ),
    frame = FALSE                  # 取消地图外边框
  ) +
  
  # # 添加地市英文名称
  # tm_text("nameE", size = 1, 
  #         root = TRUE, # 根据地图几何中心放置
  #         just = "center") +
  

  # 添加地市英文名称
  tm_text("nameE", 
          size = 1,  # 减小字体大小
          auto.placement = TRUE,
          col = "black",
          # root = TRUE, # 根据地图几何中心放置
          just = "center"
          # text.color = "black" ,
          # remove.overlap = TRUE
          ) +  # 避免文本重叠
  
  # 添加简单的指北针
  tm_compass(
    type = "arrow",  # 简单箭头指北针
    size = 2,        # 指北针大小
    position = c("right", "top"), # 指北针位置（右上方）
    text.color = "black"   # 文字颜色
  ) +
  
  # 添加比例尺
  tm_scale_bar(
    breaks = c(0, 200, 400, 600, 800), # 比例尺分段
    text.size = 0.8,         # 文字大小
    color.dark = "black",    # 深色部分颜色
    position = c("right", "bottom") # 比例尺位置（右下方）
  ) #+
  # 
  # ## 数据来源标注
  # tm_credits(
  #   text = "数据来源：中国疾病预防控制信息系统",  # 标注内容
  #   position = c("left", "bottom"),  # 左下角位置
  #   size = 0.7                      # 文字字号
  # )


## 输出交互预览
print(epidemic_map)

## 保存出版级矢量图
tmap_save(
  tm = epidemic_map,
  filename = "./主图/地市累计发病率地图.jpg",
  # width = 16,             # 图像宽度（厘米）
  # height = 12,            # 图像高度
  # units = "cm",           # 尺寸单位
  dpi = 600               # 打印分辨率
)

# 4.2 可视化设计（发病数） -----------------------------------------------------------------

## 计算2005-2023年累计发病数
cumulative_cases <- rate_year_city %>%
  filter(between(year, 2005, 2023)) %>%  # 筛选目标年份
  group_by(地区, 地区编码) %>%          # 按地州分组
  summarise(
    total_cases = sum(count, na.rm = TRUE),  # 总发病数
    .groups = "drop"
  )

## 匹配地理数据与发病数数据
geo_cases <- geo_data %>% 
  left_join(cumulative_cases, by = c("name" = "地区")) %>%  # 按名称匹配
  mutate(
    total_cases = replace_na(total_cases, 0)  # 缺失值填充为0
  ) %>% left_join(nameCE_city, by = c("地区编码" = "city_code"))

geo_cases$total_cases %>% unique()
geo_cases$total_cases %>% summary()

## 构建基础地图框架
epidemic_map <- tm_shape(geo_cases, bbox = st_bbox(geo_cases)) +  # 设置绘图范围
  
  ## 核心制图元素：发病数分级填色
  tm_polygons(
    col = "total_cases",        # 映射到发病数字段
    title = "Cumulative Cases\n(2005-2023)",
    palette = "YlOrRd",  # 使用WHO推荐的色阶
    style = "fixed",     # 使用固定间隔分级
    breaks = c(0, 1000, 3000, 6000, 9000, 12000,15000), # 手动设置分组
    border.col = "gray40", # 边界颜色
    lwd = 0.6,           # 边界线宽
    alpha = 0.8          # 透明度
  ) +
  
  ## 地图布局设计
  tm_layout(
    legend.position =  c("left", "top"),
    legend.format = list(
      fun = function(x) format(x, digits = 2)  # 图例数值保留2位小数
    ),
    frame = FALSE                  # 取消地图外边框
  ) +
  
  # 添加地市英文名称
  tm_text("nameE",  # 确保使用正确的列名
          size = 1,  # 减小字体大小
          auto.placement = TRUE,
          col = "black",
          just = "center") +  # 避免文本重叠
  
  # 添加简单的指北针
  tm_compass(
    type = "arrow",  # 简单箭头指北针
    size = 2,        # 指北针大小
    position = c("right", "top"), # 指北针位置（右上方）
    text.color = "black"   # 文字颜色
  ) +
  
  # 添加比例尺
  tm_scale_bar(
    breaks = c(0, 200, 400, 600, 800), # 比例尺分段
    text.size = 0.8,         # 文字大小
    color.dark = "black",    # 深色部分颜色
    position = c("right", "bottom") # 比例尺位置（右下方）
  )


## 输出交互预览
print(epidemic_map)

## 保存出版级矢量图
tmap_save(
  tm = epidemic_map,
  filename = "./主图/地市累计发病数地图.jpg",
  dpi = 600               # 打印分辨率
)

